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The Neel temperature, Tn, of quasi-one- and quasi- two-dimensional antiferromagnetic Heisenberg 
models on a cubic lattice is calculated by Monte Carlo simulations as a function of inter-chain (inter- 
layer) to intra-chain (intra-layer) coupling J'/J down to J'/J ~ 10^'^. We find that Tn obeys a 
modified random-phase approximation-like relation for small J' j J with an effective universal renor- 
malized coordination number, independent of the size of the spin. Empirical formulae describing 
Tn for a wide range of J' and useful for the analysis of experimental measurements are presented. 



While genuinely one-dimensional (ID) and two-dimen- 
sional (2D) antiferromagnetic Heisenberg (AFH) models 
cannot display long-range order (LRO) except at zero 
temperature [l| , weak inter-chain or inter-layer couplings, 
J', which always exist in real materials, lead to a fi- 
nite Neel temperature Tn- So far, the J'-dcpcndcncc 
of Tn was calculated by exactly treating effects of the 
strong interaction J in the ID or 2D system, but using 
mean-field approximations for the inter-chain and inter- 
layer coupling J' Recently, more advanced theories 
of the latter effects have been proposed for quasi-lD 
(QID) d 01 and quasi-2D (Q2D) %^ systems, and the 
results have been compared with the experimental ob- 
servations on QID antiferromagnets, e.g., SroCuOs 0, 
and Q2D antiferromagnets, e.g., La2Cu04 |3- In view 
of the importance of experimentally well-studied Q2D 
antiferromagnets as undoped parent compounds of the 
high-temperature superconductors, accurate and unbi- 
ased numerical results for QID and Q2D AFH models 
are strongly desired. In a recent work along this line, 
Sengupta e< al. Q have demonstrated peculiar temper- 
ature dependences of the specific heat in the quantum 
Q2D AFH model. 

Here we calculate the Neel temperature Tn as a func- 
tion of J' in fully three-dimensional (3D) classical and 
quantum Monte Carlo (MC) simulations of coupled- 
chains and coupled-layers. Our MC results on the quan- 
tum spin-S* and classical S* = oo AFH models are ana- 
lyzed by a modified random-phase approximation (RPA) 
with a rcnormalized coordination number defined by 



C(^') 



1 



j'xs(rN(J')) 



(1) 



where Xs(?^) is the staggered susceptibility of the ID or 
2D model at temperature T. 

In a simple RPA calculation Q , this quantity is just the 
coordination number Zd in the inter-chain or inter-layer 
directions: zi = 4 and Z2 = 2 for the QID and Q2D 



systems, respectively. Our main result is that C('^') eval- 
uated by Eq. with our numerically obtained Tn(J') 
and Xs(r) becomes constant 



C( J') ~ Cd = ^dZd 



(2) 



for J' < ~ O.IJ, with the constants k\ = 0.695 and 
fc2 = 0.65. These constants fc^ differ from the simple 
RPA result /c^ = 1, but the value of k\ is consistent with 
the modified self-consistent RPA theory for the quan- 
tum QID (q-QlD) model of Irkhin and Katanin (IK) 
Furthermore we find, that, within our numerical accu- 
racy, the value of kd is the same for the S = 1/2, S ~ \, 
S ~ 3/2 and S = oo, and we conjecture that kd is uni- 
versal and independent of the spin S for small J'/ J. 

We also propose empirical formulae for Tn(J') for all 
values of J' examined in the present work up to J' = J 
where corrections to the modified RPA are significant 
quantitatively. These formulae are useful in analyz- 
ing experimental results on infinite-layer antiferromag- 
nets such as Can^Sro.i5Cu02 [3i (5CAP)2CuBr4 and 
(5MAP)2CuBr4 |lO|, where they allow to determine the 
strength of the inter-chain or inter-layer coupling J' from 
experimental measurements of Tn. 

Model and numerical methods. — The Hamiltonian 
of the QID and Q2D AFH models is defined on an 
anisotropic simple cubic lattice: 

+ Jz^i,j,k • ^i,j,k+l) , (3) 

where the summation ^ ^ runs over all the lattice sites 
on a.n Lx X 1 



-^i,j,k 

Lz cubic lattice and Si,j^k is the spin 



operator at site (i, j, k). We put Jx = Jy = J' and = J 
for the QID model and = Jy = J and Jz = J' for the 
Q2D model with J > and < J' < J. For comparison, 
we also examine the classical limit S* = oo of Eq. JSJ. 
Note that JS{S -\- 1) sets the energy scale in the classical 
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TABLE I: Neel temperatures of the S = 1/2 q-QlD, c-QlD, 
5=1/2 q-Q2D, and C-Q2D AFH models, normalized by 
JS{S + 1). The result for the classical system with J' /J = 1 
is taken from Ref. H^ . 



J'/J 


T^/JS{S + 1) 




q-QlD 


c-QlD 


q-Q2D 


C-Q2D 


1 


1.2589(1) 


1.4429(1) 


1.2589(1) 


1.4429(1) 


0.5 


0.78997(8) 


0.9317(1) 


1.0050(4) 


1.1733(1) 


0.1 


0.22555(3) 


0.39551(8) 


0.6553(4) 


0.8526(1) 


0.05 


0.12171(5) 


0.28377(4) 


0.5689(2) 


0.7797(1) 


0.02 


0.05258(1) 


0.18361(3) 


0.48463(8) 


0.7115(1) 


0.01 


0.02768(1) 


0.13157(2) 


0.43515(6) 


0.6731(2) 


0.005 




0.09393(2) 


0.39513(4) 


0.6419(2) 


0.001 




0.042547(6) 


0.32571(8) 


0.5858(4) 



model. In practice, we simulate a system of unit vectors, 
which is equivalent to fixing JS{S + 1) to unity. 

The MC simulations have been performed using the 
continuous-imaginary-time loop algorithm for the 
quantum model (QMC) and the Wolff cluster algo- 
rithm [III for the classical model (CMC). The AF cor- 
relation length in each of the directions a (= x, y, z) 
are evaluated by the second-moment method |l3l | on lat- 
tices whose aspect ratio is chosen such that ^a/La does 
not depend on a in the vicinity of Tn. Explicitly, for the 
5=1/2 q-QlD systems with J'/J =0.01, 0.05, 0.1, and 
0.5, we set L,JL^ = 36 {L,, < 504), 12 (384), 4 (200), 
and 2 (128), respectively, while for the S" = 1/2 quantum 
Q2D (q-Q2D) systems with J'/J = 0.001, 0.005, 0.01, 
and > 0.02, we set L^/L, = 48 (L^ < 288), 20 (240), 16 
(192), and 1 (80), respectively. Then we determine Tn 
from finite-size scaling, looking for the best data collapse 
of£,a/La plotted versus {T-Tn)L^/'' for different system 
sizes. Wc have fixed the exponent v = 0.71 |Q| to the 
value of the 3D classical Heisenberg universality class. 
The values of Tn obtained for the S" = 1/2 q-QlD and 
classical QID (c-QlD) systems, and the = 1/2 q-Q2D 
and classical Q2D (c-Q2D) systems are listed in Tabled 

QID systems. — The classical ID (c-lD) model shows 
LRO at T = 0, while the ground state of the 5=1/2 and 
S = 3/2 quantum ID (q-lD) model is gapless and has no 
LRO [13 ■ Correspondingly, the staggered susceptibility 
Xs for the classical model, given exactly by 0| 

^ a: 1 + 1/tanh.x-l/a; 
^ 3Jl-l/tanh.T + l/x- ^ ' 

with X = JS{S+1)/T, diverges as T''^ in the limit T 
0, while the one for the 5 = 1/2 model, asymptotically 
given by 0| 

Xs{T) ^ fyin(^) + ilnln(^) , (5) 

exhibits only a 1/T divergence with logarithmic correc- 
tions. Here we note that the quantitative accuracy of this 
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FIG. 1: J'-dependence of rN/JS(S'+ 1) for the QID systems. 
The open symbols denote our MC results. The error bar of 
each point is much smaller than the symbol size. The dashed 
curve passing through the S = 00 data is obtained by Eq. Q 
with fci=0.7 and Eq. Q for Xs. The solid curve denotes the 
proposed empirical formula of Eq. The other curves are 
results of various approximations discussed in the text. 

expression is limited to a very low temperature range. In 
fact Eq. (O with the constants ci and A derived field- 
theoretically ^3 does not fit well to Xs calculated nu- 
merically at T > 0.003J. This indicates the limits of 
applicability of analytical results and that one has to use 
instead numerical data in this temperature range. 

Due to the different functional forms of the quantum 
and classical susceptibilities, we observe in Fig. ^ that, 
at small J'/J, Tn(J') oc ^ J'/J for the classical model, 
while Tj^{J') oc J'/J with logarithmic corrections for the 
quantum model. Comparing the RPA result (Eq. Q 
with fci = 1) with the modified RPA one (Eq. with 
ki ~ 0.70, denoted by IK), one can easily see that the lat- 
ter describes Tj<!{J') much better and is a fairly good de- 
scription of Tn(J') in the range J'/J < 0.3. Comparing 
our results to the next leading order finite-temperature 
perturbation theory 0| (NLO in Fig. ^ which is based 
on Eq. (jSJ, however, we do not find good agreement, be- 
cause, as pointed out above, Eq. (0) is inappropriate in 
the considered temperature range. 

The agreement with the modified RPA theory is di- 
rectly shown in Fig.|21where the J'-dependence of C( J') in 
Eq. (P) is shown. The Xs(In( J')) are obtained from QMC 
simulations interpolated near T = Tn for the 5=1/2 
model and from Eq. I^J for the 5 = 00 model. For 
J'/J < 0.1 we reach Eq. Q with ki ~ 0.695 for the 
5=1/2 model as well as for the classical limit 5 = 00 
model and conclude that, within the numerical accuracy 
of our simulation, the modified RPA with J'-independent 
C(J') is an appropriate quantitative description of the 
models with J'/J in this range. 

Interestingly the result mentioned above seems to hold 
for quantum models with other values of 5. As also 
shown in Fig. [3 within our numerical accuracy, this is 



3 





0.74 




0.72 




0.7 








0.68 




0.66 




0.64 




0.62 



QID 

-A- -A A- 

5=1/2 (QMC) 
1 

3/2 
oo (CMC) 



□ 

O [] 



0.001 



0.01 



0.1 



J' IJ 



FIG. 2: J'-dependence of C^{J')/z\ for the QID systems. In 
all cases ({J')/zi approaches a constant (~ 0.695), denoted 
by the dotted line, at small J' /J. The error bar of each point 
is smaller than the symbol size. 



well confirmed for the S = 3/2 model with J' /J > 0.02. 
For the 5=1 case we find agreement in the range 
jyj > 0.05, where Tjv is larger than the Haldane 
gap 01 of the isolated chain. Below this temperature 
the finite size scaling of the QMC data becomes less reli- 
able and we cannot draw definitive conclusions. Even if 
the result for the S' = 1 model is restricted to this tem- 
perature range, the present result is surprising, given the 
different behavior of Xs{T) in the c-lD and q-lD models. 

Q2D systems. — In both classical 2D (c-2D) and quan- 
tum 2D (q-2D) models, AF-LRO appears at T = 0, to- 
gether with an exponential divergence of XsiT) at T ^ 0. 
In the C-2D system, Ys is proportional to exp(47r J/T) 
at low temperatures [ijj and our numerical results agree 
with previous simulations For the q-2D models, 

there is a similar exponential divergence at T ^ 0. 
In the renormalized classical regime of the non-linear a 
model ji^l, for example, Xs{T) is written as 



Xs(T) J = C2T/J exp(47rp,/T), 



(6) 



where ps is the spin stiffness and C2 a constant. 

The J'-dependence of Tn for the Q2D models is shown 
in Fig. El We see that Tj^{J') oc -I/ln(J7J) at small 
J'/ J in the S = 1/2, 5 = 1 and S = oo models due 
to the similar exponential forms of Xs at T of the 
classical and quantum models. Figure0]shows that again 
for J' /J < 0.05 the values of C{J') are universal for the 
quantum and the classical models: k2 = 0.65 in Eq. (O 
independent of the spin size S. This confirms the validity 
of our modified RPA scenario represented by Eqs. and 
(121) also for the Q2D systems. 

If we insert Eq. ^ into Eq. (P) with ({J') = C2, we 
obtain the following expression of Tn for J' /J <C 1, 

Tn = 4^p,/(& - In ( J'/J) - In (Tn/ J)) , (7) 

with b = — ln(C2C2). This result is compatible with that 
of the 1/A^-cxpansion theory of Tn(J') due to Irkhin et 




J'/J 

FIG. 3: J'-dependence of Tn/JS{S + 1) for the Q2D sys- 
tems. The open symbols denote our numerical results for the 
S = 1/2 (circles), S — 1 (squares), and S — 00 (triangles) 
models. The error bar of each point is much smaller than the 
symbol size. The dashed (RPA) and dotted (modified RPA) 
curves for 5 = 1/2 are obtained from Eqs. and (|2j with 
k2 — 1 and k2 = 0.65, respectively. The solid curves denote 
the proposed empirical formula @, while the curve passing 
through the 5 = 00 data is simply a guide for the eye. The 
inset shows the same data on a linear scale. 



al. 13 for the 5=1/2 model in the same limit. Various 
estimations of b and ps can be obtained analytically 
according to the different approximation schemes used. 
Unfortunately, we cannot judge which approximation is 
most relevant in general since higher order corrections in 
T/Airps over the leading asymptotic expression Eq. Q 
are known to be necessary |22l | to reproduce the numeri- 
cally obtained Xs in the temperature range T /A-nps > 0.1 
simulated. In fact, corrections of this type and uncer- 
tainty on T/v due to the different possible estimates of b 
are comparable. We expect, on the other hand, that the 
constancy of the normalization factor k2 , which is found 
numerically to be within 2% in 0.001 < J'/J < 0.05 and 
0.32 < Tn/J5(5+ 1) < 0.57 (Fig. 4), holds in the limit 
J'/J ^ as well. 

Empirical formulae. — Finally, we propose empirical 
formulae for Tj^{J') based on our QMC results. For the 
5=1/2 q-QlD system we propose a modified RPA form 
based on Eqs. Q and jSJ with a constant 



j' = TN/(4c\/ln(^) + ilnln(^)), 



(8) 



but with modified values of the constants with c = 0.233 
and A = 2.6. These values are chosen to reproduce not 
Xs(T) but Tn(J') in a wide range of J'. This formula 
describes Tn(J') very well in the whole range of J'/J 
as shown in Fig. ^ and it can be used to analyze experi- 
mental results, e.g. to obtain J'/J ~ 0.0007 for Sr2Cu03 
from Tn/J~ 0.002 0. 

For the q-Q2D systems, we find that instead of Eq. {TJ, 
the following simpler expression describes Tn better in 
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FIG. 4: J'-dependence of C,{J')/z2 for the Q2D systems. In 
all cases ({J')/z2 approaches a constant (~ 0.65), denoted by 
the dotted line, at small J'/ J. The error bar of each point is 
smaller than the symbol size unless given explicitly. 

TABLE II: Inter-layer coupling J' estimated by Eq. Q for 
various infinite- layer compounds. The Neel temperatures Tn, 
the intra-layer couplings J, and their ratio estimated by the 
experiments are also listed. 



Compound 


Tn 


J 


Tn/J 


J'/J 


Cao.ssSro.isCuOa [2] 


537K 


1535K 


0.35 


0.016 


(5CAP)2CuBr4 [10] 


5.08K 


8.5K 


0.60 


0.24 


(5MAP)2CuBr4 ]10} 


3.8K 


6.5K 


0.58 


0.22 


(5CAP)2CuCl4 [10] 


0.74K 


1.14K 


0.64 


0.33 


(5MAP)2CuCl4 [10] 


0.44K 


0.76K 


0.57 


0.21 



the range 0.001 < J'/J < 1 (see Fig.^: 
Tn = 47rp,/(fo- In (J'/J)) 



(9) 



with ps/J = 0.183 and b = 2 A3 for S = 1/2, and p^/J = 
0.68 and 6 = 3.12 for S" = 1. Table IHl shows inter-layer 
couplings J' estimated using this equation for a number 
of infinite-layer materials with S = 1/2. 

To conclude, we have determined, by high-precision 
Monte Carlo simulations, the Neel temperatures of quan- 
tum and classical QID and Q2D Heisenberg antiferro- 
magnets. Besides finding empirical formulae for Tn(J'), 
we observe that, using numerically accurate values of 
the staggered susceptibility a modified RPA with the J'- 
independent renormalized coordination number (^d suc- 
ceeds in quantitatively describing the relation between 
Tn and J'/J for J' < J^ with J^ ~ 0.1 J. 

An intriguing result of our simulations is the indepen- 
dence of on the value of the spin S, suggesting a uni- 
versality of corrections to RPA for J' ^ J. Since in this 
temperature regime the physics of all these models should 
be well described by an anisotropic non-linear cr-model 
(NLctM) in the renormalized classical regime, we conjec- 
ture universal corrections to RPA also for the NLcrM. 
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